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Abstract 

We present a model of solids made from polygonal cells connected 
via beams. We calculate the macroscopic elastic moduli from the 
beam and cell parameters. This modellisation is particularly suited 
for the simulation of fragmentation processes. We study the effects of 
an explosion inside a circular disk and the impact of a projectile and 
obtain the fragment size distribution. We find that if breaking only 
happens under tensile forces a layer on the free wall opposed to impact 
is first ejected. In that case the distribution follows a power-law with 
an exponent that in most cases is around two. 



1 Introduction 

Fragmentation plays an important role in a wide variety of physical phe- 
nomena. Examples range from geophysics to astrophysics: fragments from 
weathering, coal heaps, rock fragments from chemical and nuclear explo- 
sions, projectile collisions, asteroids etc. Most of the measured fragment 
size distributions show power law behavior with exponents between 1.9 and 
2.6 concentrating around 2.4 - 0|. Power law behavior for small fragment 
masses seems to be a common characteristic of the breaking of brittle objects. 

Fragmentation in one dimension was studied using thin glass rods dropp- 
ed vertically onto the floor [Q. Depending on the height from which the rod 
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was dropped the fragment size distribution varies from log - normal to power 
law with increasing height. 

Recently, Oddershede, Dimon and Bohr reported that the mass distri- 
bution of fragments produced by impact experiments shows universal power 
law behavior. The scaling exponents depend on the overall morphology of 
the objects but are independent of the type of the material. (Their experi- 
ment was performed with several different kinds of materials, gypsum, soap, 
paraffin, potato.) Because of the universal power law behavior found without 
any control parameter, they concluded that the fragmentation process is a 
self - organized critical process with exponent ~ 1.6 

Several theoretical approaches were established to describe fragmentation. 
In one dimensional stochastic models power law, exponential and log-normal 
fragment mass distributions can be obtained depending on the fracture point 
distributions 

Discrete stochastic processes have also been studied as models for frag- 
mentation using cellular automata. In Ref. 8 two- and three - dimensional 
cellular automata were proposed to model power law distribution in shear 
experiments on a layer of uniformly sized fragments. 

In Ref. 9 an iterative stochastic process is studied as a model of two- and 
three - dimensional discrete fragmentation. Depending on the parameters 
of the model, log - normal and power law distribution were found for the 
fragment size distribution. 

The mean - field approach describes the time evolution of the concentra- 
tion c{x, t) of fragments having mass x through a linear integro - differential 
equation |]10| : 

^ ' ' = -a{x)c{x,t) + c{y,t)a{y)f{x\y)dy (1) 



dt 

where a{x) is the overall rate at which x breaks in a time interval dt, while 
f{x\y) is the relative rate at which x is produced from the break-up of y. 
With some further assumptions on f{x\y) exact results can be obtained but 



in physically interesting situations the solution is very difficult [10 



Three - dimensional impact fracture processes of random materials were 



modeled based on competitive growth of cracks |]Tl|. A universal power 
law fragment mass distribution was found consistent with self - organized 
criticality with an exponent of |. 

Recently, a two dimensional dynamic simulation of solid fracture was 



performed using a cellular model material [T^, 0. The compressive failure 



of a rectangular sample, the four - point shear failure of a beam and the 
impact of particles with a plate and with other particles were studied. 
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In this paper we present a two dimensional dynamic simulation of the 
fragmentation of granular materials. We establish a two - dimensional model 
for a deformable, breakable, granular solid by connecting unbreakable, un- 
deformable elements with elastic beams, similar as in Refs. 12,13 . The 
contacts between the particles can be broken according to a certain breaking 
rule, which takes into account the stretching and bending of the connections. 
The breaking rule contains two parameters to describe the relative impor- 
tance of the two breaking modes. We address the question how the fragment 
mass distribution depends on the breaking parameters. After measuring the 
macroscopic elastic behavior of our granular solid the model is apphed to 
study fragmentation caused by an explosion inside a solid and by a projectile 
impact. 



2 The model 

In order to study the fragmentation of granular solids we use Molecular 
Dynamics (MD) simulations in two dimensions. A general overview of MD 
simulations applied to the field of granular materials can be found in Ref. 14. 
This method calculates the motion of particles by solving Newton's equations. 
In our simulation this is done using a Predictor Corrector scheme. 

Our model of a deformable, breakable granular solid is a generalization 
of the model which was used earlier to study the flow of granular materials 
1l5| . The model construction is composed of three major steps, namely, the 



implementation of the granular structure of the solid, the determination of 
the elastic behavior by the contact force and the beam models, and finally 
the breaking of the solid. In this section, the three steps of the model con- 
struction are presented in detail. 



(a) Granularity 

In order to take into account the complex structure of the granular solid 
we use arbitrarily shaped convex polygons. To get an initial configuration of 
these polygons we construct a vectorizable random lattice, which is a Voronoi 
construction with a regularization procedure (see Ref. 16). The advantage of 
the vectorizable random lattice compared to the ordinary Poissonian Voronoi 
tessellation is that the number of neighbors of each polygon is limited which 
makes the computer code faster and allows to simulate larger systems. The 
convex polygons of this Voronoi construction are supposed to model the 
grains of the material, see Ref. 15. In this way the structure of the solid is 
built on a mesoscopic scale. Each element is thought of as a large collec- 
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tion of atoms. In our simulation however these polygons are the smallest 
particles interacting elastically with each other. All the polygons have three 
continuous degrees of freedom in two dimensions: the two coordinates of the 
positions of the center of mass and the rotation angle. 

(b) Elastic behavior of the solid 



The elastic behavior of the solid is captured in the following way: The 
polygons are considered to be rigid bodies. They are not breakable and not 
deformable. But they can overlap when they are pressed against each other 
representing to some extent the local deformation of the grains. Usually the 
overlapping polygons have two intersection points. These points define the 
contact line (see Fig. |I]). 

In order to simulate the elastic contact force between touching grains we 
introduce a repulsive force between the overlapping polygons. This force is 
proportional to the overlapping area A divided by a characteristic length 
Lc of the interacting polygon pair. Our choice of Lc is given by 1/Lc = 
l/2(l/rj + l/vj), where ri,rj are the diameters of circles of the same area 
as the polygons. This normalization is necessary in order to reflect the fact 
that the spring constant is proportional to the elastic modulus divided by a 
characteristic length. (In the case of a linear spring this characteristic length 
is simply the equilibrium length of the spring.) The direction of the force is 
chosen to be perpendicular to the contact line of the polygons. 

The contact force Fij between two particles is given by 



YA 

= -—n (2) 



c 



where n is the unit vector perpendicular to the contact line (see Fig. [l|) and 
Y is the grain bulk Young modulus. The friction of the touching polygons 
can be implemented according to Coulomb's friction law (see also Ref. 15). 
However it turned out from the simulations that friction does not play an 
important role in the fragmentation of the granular solids. Therefore, in the 
present simulations the friction term was omitted. 

In order to keep the solid together it is necessary to introduce a cohesion 
force between neighboring polygons. For this purpose we introduce beams, 
which were extensively used recently in crack growth models |T^, 113. The 



centers of mass of neighboring polygons are connected by elastic beams, 
which exert an attractive, restoring force between the grains, and can break 
in order to model the fragmentation of the solid. 
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Because of the randomness contained in the Voronoi - tessellation the 
lattice of beams is also random. An example of a random lattice of beams 
coupled to the Voronoi polygons can be seen in Fig. |^. 

A beam between sites i and j is thought of as having a certain cross 
section S^^ giving to it not only longitudinal but also shear elasticity. This 
cross section is the length of the common side of the neighboring polygons in 
the initial configuration. The length of the beam is defined by the distance 
of the centers of mass. The elastic behavior of the beams is governed by two 
material dependent constants. For the beam between sites i and j: 

a-' = ^ (3) 

b'^ = (4) 
GS'3 ^ ' 



EFi 



(5) 



where E and G are the Young and shear moduli of the beam, S^^ is the 
area of the beam section, and is the moment of inertia of the beam for 
flexion. A flxed value of E was used for all the beams and was chosen to 
be = 2a^K The length, the cross section and the moment of inertia of each 
beam are determined by the random initial conflguration of the polygons as 
explained above. The beam Young modulus E and the grain bulk Young 
modulus Y are, in principal, independent. 

In the local frame of the beam three continuous degrees of freedom are 
assigned to both lattice sites (centers of mass) connected by the beam, which 
are for site i, the two components of the displacement vector and a 

bending angle 6*. 

For the beam between sites i and j one has the longitudinal force acting 
at site i: 



the shear force 



(6) 



f; = pi^(ui-ui)-t^{Q^ + en, (7) 



and the flexural torque at site i 

Mi = ^^{ui~ui + pe^) + s'n'^\e^ -Q') (s) 

2 

where a'^ = p'^ = l/{b'^ + l/12d^), and 6'^ = p'^{b'^ /d^ + 1/3) . It can 

be shown that this beam model is a discretisation of the simplified Cosserat - 
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equations of continuum elasticity which should be used to describe the elastic 
behavior of the granular solids, instead of the Lame equations [l^ . 



( c ) Breaking of the solid 

To model fragmentation it is necessary to complete the model with a 
breaking rule, according to which the over-stressed beams break. 

For not too fast deformations the breaking of a beam is only caused by 
stretching and bending. We impose a breaking rule which takes into account 
these two breaking modes, and which can reflect the fact that the longer and 
thinner beams are easier to break. We used a breaking rule of the form of 



the von Mises plasticity criterion [T^: 



2 



max{\<d'^\, 



+ -^^7 '"^ > 1 (9) 

where e = Al/l is the longitudinal strain of the beam, and 0^ are the 
rotation angles at the two ends of the beam and and are threshold values 
for the two breaking modes. In the simulations we used the same threshold 
values and te for all the beams. 

Two sets of simulations were performed applying breaking criteria of the 
type of Eq. (|^) but once only for stretched beams, i.e. e > and once for 
stretched and compressed beams , i.e. for all e. The first case is physically 
more relevant since it reflects the fact that it is much harder to break a solid 
under compression than under elongation. 

The first term of Eq. takes into account the role of stretching and the 
second term the role of bending. Varying the threshold values the relative 
importance of the two modes in the beam breaking can be changed. 

During the simulation the left hand side of Eq. (y) is evaluated at each it- 
eration time step for all the existing beams, which fulfill the strain conditions. 
The breaking of beams means that those beams for which the condition of 
Eq. holds are removed from the calculation, i.e. their elastic constants 
are set to zero. Removed beams are never restored during the simulation. 

The surface of the grains, on which beams are broken represents cracks. 
The energy of the broken beams is released in creating these new crack sur- 
faces inside the solid. 



3 Results 

With the model described above it is possible to perform a variety of exper- 
iments. 
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First of all it is necessary to study the global elastic behavior of our 
two dimensional solid assembled from cells. Simulations were performed to 
measure the Young modulus and the Poisson ratio of a rectangular sample. 

In the present paper we are mainly interested in fragmentation. We use 
the model to examine the fragmentation of a heterogeneous solid in two 
different experimental situations, namely, the fragmentation of a disc-shaped 
solid caused by an explosion, and the fragmentation of a rectangular solid 
block due to the impact of a high velocity projectile. 

3.1 Elastic properties 

In Ref. 12 it was shown that the macroscopic elastic behavior of two dimen- 
sional materials composed of particles with linear contacts can be character- 
ized by two elastic constants, the Young modulus K and the Poisson ratio 
similar to the case of homogeneous isotropic solids. It seems clear from the 
definition of the model that these elastic constants depend on the properties 
of the constituent particles, i.e. on the shape of the grains, on the stiffness 
of the grain - grain contacts (in our context the grain bulk Young modulus 
Y and the beam Young modulus E) and on the typical grain size. 

For real materials the Young modulus K and the Poisson ratio v are 
usually determined by uniaxial loading of the body. In order to measure 
numerically the elastic properties of our simulated material we apply uniaxial 
loading on a two dimensional rectangular sample, which has linear extensions 
L in the direction of the loading and S in the perpendicular direction. The 
corresponding changes of the extensions due to loading are denoted by AL 
and /S.S. During the simulations plane - stress conditions are used, i.e. it 
is assumed that there is no stress in the direction out of the plane of the 
material. The boundaries parallel to the loading are free in the simulation. 

In the case of uniaxial loading of a rectangular sample under plane - stress 
conditions, the Young modulus K is the ratio between the external force F 
per unit length acting on the solid and the strain in the direction of the 
loading, F/ S — K /S.L/ L. The Poisson ratio v is the ratio between the strain 
in the direction perpendicular to the loading and the strain in the direction 
of the loading, AS/S = —uAL/L. 

To avoid the disturbing effect of the elastic waves induced by the loading, 
the numerical experiment is performed in the following way (see also Ref. 12): 
The two opposite boundaries of the solid start to move with zero initial 
velocity and non-zero acceleration. When a certain velocity is reached the 
acceleration is set to zero and the velocity of the boundaries is kept fixed. 
With this slow loading the vibrations of the solid can be reduced drastically 
compared to the case when the boundaries start to move with non - zero 
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initial velocity. A further way to suppress artificial vibrations is to introduce 
a small dissipation (friction or damping) between the grains. This dissipation 
has to be small enough not to affect the quasi-static results. To obtain the 
values of K and the force acting on the boundary layer, and the horizontal 
and vertical extensions of the sample were monitored during the loading. 

Fig. shows an example of the results of the measurement of the Young 
modulus, the horizontal stress a = F/S as a function of the vertical strain 
e = AL/L. The values of K were extracted from the slope of the straight 
lines. It can be seen that elongation and compression of the system are not 
symmetric. The solid is more stiff under compression. The reason of the 
asymmetry is that in the case of elongation practically only the beams act 
but under compression one measures the common effect of the beams and 
the overlap force giving rise to a larger effective Young modulus. 

Calculations were performed for elongation and compression of the sample 
fixing the value of the grain bulk Young modulus Y and varying the beam 
Young modulus E. The results for K are shown in Fig. ^. In the case 
of elongation K is a linear function of E as expected. For large E the 
elongation and compression curves are parallel. The difference between them 
is determined by the grain bulk Young modulus Y and by the geometry of 
the Voronoi tessellation. Small values of E compared to Y means that the 
cohesive force in the solid is small with respect to the repulsive overlap force. 
Under compression the polygons can move perpendicular to the loading in the 
direction of the free boundaries in the limiting case of small cohesion. Because 
of the dense packing (the initial packing fraction is one) the deformation tends 
to increase the overall volume resulting in dilatancy [^. This effect gives 
rise to a small effective Young modulus, such that K is even smaller than 
the grain bulk Young modulus Y. 

For the Poisson ratio u a similar asymmetry of elongation and compression 
is observed as shown in Fig. ^. The values of u for elongation are smaller 
than for compression but for large E both values approach the Poisson ratio 
of the pure beam lattice. In the limiting case of small cohesion i> can even 
exceed unity, due to the dilatancy mentioned above. 

3.2 Explosion of a disc - shaped solid 

The catastrophic fragmentation of solids will be studied in two different ex- 
perimental situations: through an explosion which takes place inside the 
solid and through the impact with a projectile (stroke with a hammer). This 
chapter is devoted to study the explosion. 

In the explosion experiment the detonation takes place in the center of 
a solid disc. The granular solid with disc - like shape was obtained starting 
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from the Voronoi tessellation of a square and cutting out a circular disc in 
the center, see Fig. 2. 

In the center of the solid we choose one polygon, which plays the role 
of the explosive. Initial velocities are given to the neighboring polygons 
perpendicular to their common sides with the central one. The sum of the 
initial linear momenta has to be zero, reflecting the spherical symmetry of 
the explosion. From these two constraints it follows that for a polygon having 
mass m and a common boundary of length S with the explosive center, the 
initial velocity is proportional to ^. The sum of the initial kinetic energies 
defines the energy Eo of the explosion. (For the parameter values and the 
initial conditions of the simulation see Table |l|.) 

As a result of these initial conditions a circularly symmetric outgoing 
compression wave is generated in the solid. In our context this means that 
there is a well - defined shell where the average longitudinal strain of the 
beams < e >=< Al/l > is negative. This compression wave is not homo- 
geneous in the sense that not all the beams in this region are compressed . 
If the angle of a beam with respect to the radial direction is close to 7r/2 a 
beam can be slightly elongated within the compression wave. 

Since the overall shape of the solid has the same symmetry as the com- 
pression wave it is possible to avoid geometrical asymmetries, which would 
arise for example in the explosion of a rectangular sample due to the corners. 
In the disordered solid the initial compression wave gives rise to a compli- 
cated stress distribution, in which the over-stressed beams break according 
to the breaking rule Eq. (^). The simulation is stopped if there is no beam 
breaking during 300 successive time steps. Free boundary conditions were 
used in all simulations. 

Due to the beam breaking the solid eventually breaks apart, i.e. at the 
end of the process it consists of well separated groups of polygons. These 
groups of polygons, connected by the remaining beams, are the fragments. In 
the simulation of the explosion we are mainly interested in the time evolution 
of the fragmentation process and the mass distribution of fragments at the 
end of the process as a function of the breaking parameters. 

All the calculations were performed on the CMS of the CNCPST in Paris. 
We used the farming method, i.e. the same program runs on a variety of 
nodes with different initial setups. In our case 32 nodes were used with 
different seeds for the Voronoi generator, i.e. with differently shaped Voronoi 
cells. 
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3.2.1 The time evolution of the explosion 

We performed two sets of calculations one allowing the beam breaking solely 
under elongation and another under elongation and compression, keeping all 
the other parameters of the simulation fixed. The evolution of the explosion 
and the resulting breaking scenarios are different in these two cases. 

We can distinguish two regimes in the time evolution of the explosion: 
The initial regime is controlled by the compression wave and the disorder of 
the solid. The amplitude of the shock wave is proportional to the ratio of 
the average initial velocity of the polygons to the longitudinal sound speed 
of the solid. The width and the speed of the wave are mainly determined by 
the grain size and the Young moduli. 

If the beams are not allowed to break under compression the compression 
wave can go outward almost unperturbed and an elongation wave is formed 
behind it. The elongated beams break according to the breaking rule Eq. 
d^) only for e > 0. Due to the elongation wave a highly damaged region is 
created in the vicinity of the explosive center where practically all the beams 
are broken and all the fragments are single grains. This highly damaged 
region is called the mirror spot. Since the breaking of the beams, i.e. the 
formation of cracks in the solid, dissipates energy after some time the growth 
of the damage stops. The size of this mirror spot is determined by the 
initial energy of the explosion, by the dissipation rate and by the breaking 
thresholds (see Fig. 1). 

During and after the formation of the mirror spot when the outgoing 
compression and elongation waves go through the solid the weakest (i.e. the 
longest and thinnest) beams break in an uncorrelated fashion creating iso- 
lated cracks in the system. The uncorrelated beam breaking is dominated by 
the quenched disorder of the solid structure. This first uncorrelated regime of 
the explosion process lasts till the compression wave reaches the free bound- 
ary of the solid. 

From the free boundary the compression wave is refiected back with oppo- 
site phase generating an incoming elongation wave. The constructive inter- 
ference of the incoming and outgoing elongation waves gives rise to a highly 
stretched zone close to the boundary. The beams having small angle with re- 
spect to the radial direction have the largest elongation. In this zone a large 
number of beams break causing usually the complete break-off of a bound- 
ary layer along the surface of the solid. The thickness of this detached layer 
is roughly half the width of the incoming elongation wave (see Fig.^). The 
fragments of this boundary layer fly away in the radial direction with a high 
velocity carrying with them a large portion of the total energy in the form 
of their kinetic energy. After that the system starts to expand. This overall 
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expansion initiates cracks going from inside to outside and from outside to 
inside. The branching of single cracks and the interaction of different cracks 
give rise to the final fragmentation of the solid (Fig.^. This second part of 
the evolution of the explosion process is dominated by the correlation of the 
cracks. 

The propagation of the elastic waves when the beam breaking is switched 
off is presented in Fig. ^ One can observe the peak of the initially imposed 
shock, the propagation of the compression and elongation waves and the 
formation of the highly stretched zones at the boundary. 

If the beams are allowed to break also under compression the time evolu- 
tion of the process is significantly different. In the vicinity of the explosion 
center the initial compression wave governs the formation of the mirror spot. 
Where the compression wave passes all the beams break as long as the am- 
plitude of the compression wave is larger than the threshold strain of the 
breaking criterion Eq. (^. Thus when the energy dissipation stops the 
growth of the damage the amplitude of the remaining compression wave is 
weaker and the size of the mirror spot is larger than it was in the previous 
case (Fig. H). A significant part of the initial energy of the explosion is dis- 
sipated by the beam breaking and is transfered to the kinetic energy of the 
single polygons. So, less fragmentation can occur in the outer regions of the 
disc as compared to the case where breaking only occurs under elongation. 

The compression wave reflects back from the free boundary but the con- 
structive interference does not cause substantial damage along the surface 
(Fig.^). That is why the final breaking picture is different from the previous 
case. 

3.2.2 The fragment mass distribution 

Fig. ^ and Fig. |^ show that the final breaking scenarios for the two break- 
ing criteria are significantly different. When the beams are also allowed to 
break under compression there are only a few small fragments (apart from 
the single polygons) and a few much larger ones. This results in a fast decay 
of the fragment size distribution. The resulting fragment mass histogram is 
presented in Fig. ^. F{m) denotes the number of fragments with mass m 
divided by the total number of fragments. As can be seen from the logarith- 
mic plot in Fig. 1^ a, an exponentially decaying function can be a reasonable 
fit to these results for approximately one order of magnitude in mass. The 
quality of the fit is demonstrated in Fig. ^b using a semi-logarithmic plot. 

The application of Eq. (|^) for e > only, is physically more relevant. 
In the following we present results of the mass distribution of the fragments 
as a function of the breaking parameters allowing the beam breaking solely 
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under stretching. Changing and t@ in the breaking criterion Eq. one 
can vary the relative contribution of the stretching and bending breaking 
modes, which can affect the crack formation and the final breaking scenario. 
The question is whether the mass distribution of the fragments is invariant 
under a change of and t^. 

We performed simulations alternatively fixing and t@ and changing the 
value of the other one, keeping all the other parameters of the simulations 
fixed. In both cases the fraction of broken beams p was monitored as a 
function of s//, where s and / are the cross section and the length of a beam, 
respectively: 

Pis/l) = b \''^;jif\ (10) 

where N{s/l) and Nhroken{s/l) denote the number of beams in the sample 
and the number of broken beams having a ratio oi s/l. 6 is a normalization 
constant. Fig. ^ and Fig. |ll] show the results for fixed te = 4° varying 



the threshold strain from 1% up to 6% and for fixed = 3% varying the 
threshold angle from 1° up to 7°, respectively. The curves are normalized 
in such a way that their integrals are equal to the overall fraction of broken 
beams po, i.e. J p{s/l)d{s/l) = Po- Hereto = Nhroken/N, where N and Nbroken 
denote the total number of beams in the sample and the total number of 
broken beams, respectively. One can observe that p{s/l) is a monotonically 
decreasing function for all parameter values, which shows that longer and 
thinner beams are always easier to break, as expected. In Fig. ^ for small 
and in Fig. ^ for small Iq the solid is more fragile, so it undergoes strong 
damage. This results in a larger po and a rather fiat breaking fraction p{s/l). 
The curves belonging to increasing and te lie below each other. Since 
by increasing a breaking threshold, the contribution of the corresponding 
breaking mode becomes less important, the fixed mode starts to dominate 
the breaking and determines the limiting curve. The results obtained by 
switching completely off one of the breaking modes confirm this argument. 



The fragment mass histograms are presented in Fig and Fig. |T3|. The 
lower cutoff of the histograms is determined by the size of the unbreakable 
polygons (smallest fragments) while the upper cutoff is given by the finite 
size of the system (largest fragment). Larger values of the overall breaking 
fraction p^ entails that the system is broken into smaller pieces. 

The histograms follow a power law for practically all the parameter pairs 
for at least one order of magnitude in mass, such that we seem to have 

F(m)^am-'^. (11) 
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The effective exponents (3 were obtained from the estimated slopes of the 
curves. The results are presented in Table |[ Apart from the case of ex- 
tremely small breaking parameters the exponent j3 only slightly varies around 
/9 = 2.0, indicating a more or less universal behavior within the accuracy, 



with which (3 was determined (±0.05). Fig. |T4| shows that increasing 
and te a crossover appears in F{m), i.e. there are two exponents. As we 
have mentioned, the initial compression wave gives rise to the break-off of 
a boundary layer, the thickness of which is determined by the wavelength 
and by the breaking parameters. For larger and te, when the system is 
tougher, a detachment of the fragmentation of the boundary layer and the 
bulk of the solid appears. The boundary gives the main contribution to the 
small fragments (containing ~ 2 - 6 polygons) and the fragmentation of the 
bulk dominates the range of the larger ones. In these cases the exponents f3 
in Table ^ were obtained in the limit of large fragments. 



3.3 Impact of a projectile 

Besides the explosion a catastrophic fragmentation of solids can also be gen- 
erated by an impact with a projectile. In nature one of the most spectacular 
examples is that the size distribution of meteorites and asteroids shows power 
law behavior. These objects are believed to have originated from the fracture 
of primitive planets due to collisions [jl], 

We applied our model to study the fragmentation of a rectangular solid 
block (e.g. a block of concrete) due to an impact allowing the beam breaking 
solely under stretching e > 0. One polygon at the lower middle part of the 
block is given a high velocity directed inside the block simulating an elastic 
collision with a projectile. The boundary conditions and the stopping condi- 
tion were the same as in the explosion experiment. The breaking thresholds 
were chosen = 3% and te = 4°. 

The evolution of the fragmenting solid block is presented in Fig. [T^. As 
in the case of the explosion, the initially generated compression wave plays 
a significant role. Since the energy of the collision is concentrated around 
the impact site of the projectile the damage is the largest in that region. 
The completely destroyed zone, where all the beams are broken stretches 
inside the solid in the forward direction resulting in the break-up of the 
solid. When the shock wave reaches the boundary at the side of the solid 
opposite to the collision point it gives rise to the break-off of a boundary 
layer. The fragments of this layer fly away in the forward direction with a 
high velocity. Some small fragments from the vicinity of the collision point 
are scattered backward. The damage in the direction perpendicular to the 
projectile is not strong, the broken boundary layer is thicker and the speed 



13 



of the fragments is smaller. 

Results of laboratory experiments on high velocity impacts can be found 
in Refs. 3,4. In Ref. 4 a picture series obtained by a high speed camera is 
presented showing the time evolution of an impact experiment. Our results 
are qualitatively in good agreement with the experimental observations. 

The resulting fragment mass histogram F{m) is presented in Fig. 
Similarly to the explosion experiment, F{m) shows power law behavior for 
approximately one order of magnitude in mass. The value of the effective 
exponent is /5 = 1.98 ± 0.05. 



4 Conclusion 

We have studied shock fragmentation in two dimensions using a cell model 
similar to the one introduced by Potapov et al |l^, |13| . We found a qualita- 



tive difference in the breaking process when breaking occurs under traction 
or under both, traction and compression. If only traction force can pro- 
duce cracks a power-law distribution of the fragment sizes is evidenced with 
a rather universal exponent around two in good agreement with most ex- 
perimental observations. Interesting is the effect of lamination in which a 
surface layer half the width of the elastic wave is detached and ejected with 
high velocity. This effect has been observed experimentally in high speed 
impacts. 

With the present study we have evidenced that the full treatment of elas- 
ticity and wave propagation are necessary if one wants to reproduce subtleties 
of fragmentation as lamination or the dependence on compressive breaking. 
Still our study makes a certain number of technical simplifications which 
might be important for a full quantitative grasp of fragmentation phenom- 
ena. Most important seems to us the restriction to two dimensions, which 
should be overcome in future investigations. The existence of elementary, 
non-breakable polygons restricts fragmentation on lower scales and hinders 



us from observing the formation of powder of a shattering transition |jTO . 

An advantage of our model with respect to most other fragmentation 
models is that we can follow the trajectory of each fragment, which is of- 
ten of big practical importance and that we know how much energy each 
fragment carries away. The polygonal structure of our solid allows us to re- 
alistically model granular or polycrystalline matter, and considering as well 
cell repulsion as beam connectivity gives us a rich spectrum of possibilities 
ranging from breaking through bending to the effect of dilatancy. If one or 
the other mechanism is turned off we have the extreme cases of an elastic 
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homogeneous solid and a compact dry granular packing. 
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Parameter 


Symbol 


Unit 


Value 


Density 


P 


g/cm^ 


5 


Grain bulk Young modulus 


Y 


dyn / cm^ 


IQio 


Beam Young modulus 


E 


dyn/cm^ 


5-109 


Time step 


dt 


s 


IQ-^ 


Diameter of the solid 


d 


cm 


40 


Average # of polygons 


n 




1100 


Energy of the explosion 


Eo 


erg 


5-109 


Average initial speed 


Vo 


m/ s 


200 


Estimated sound speed 


c 


m/ s 


900 



Table 1: The parameter values used in the simulations. 





= 3% 






te 




1 % 


(4.34) 


1° 


(4.42)'^ 


2 % 


2.39 


2° 


2.5 


3 % 


2.01 


3° 


1.98 


4% 


1.97 


4° 


2.01 


5 % 


1.95 


5° 


1.96 


6 % 


1.96 


6° 


1.95 






7° 


1.97 



Table 2: The values of the exponent (3 for different breaking thresholds. 

"These exponents belong to the hmiting case of extremely small breaking thresholds. 
They were extracted by fitting straight lines over the whole mass range of the corresponding 
curves. 
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Figure 3: 



'ypical result of a measurement of the effective Young modu" 



tJ^. |hip]tii§3 ntal streag = F/S ba a function of the k^ertifcal atlrain t: 



us K, 
AL/L. 



Positi-vt^. m li^gcfti<^gtt-aA. ^fi^elo^atioi^a&a 0:641^^1 o^^tfi^ SUlfe, 
respectively. The values of K were extracted from the slopes of the straight 
lines. One can observe the asymmftey Bl cj'igation and compression. 
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Figure 4: Thj. effective (macroscopic) Young modulus K of the grs.nular 
solicb-wtfeasmpn] under rrhmprPSHinn RTlirl fflnnigatrinn I varying the Iheam roung 



moduluslfifeeSveen HO^^^P/cm^ aidi4-Q0^°(iy^/emf OThe gdaciaMlk Young 
modulus was fixed to be y = ltf°dup/cmj. 

11 {uyu/cm. ) 
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FiguK! 5^Tbe-B&is*ni"ratio of "Che granular solid measured under compression 
aQf^Tongation iThd^ valuesiof E arid Y, are, tha saitie as iniFigi ^1 The curve 



ofltheCtflare bealnc -lattOice was 2>.lit€iiiOEd by s^aitciiini^ off the-layciUSp force 
between the polygons. n / i r 

h (dyn/cm } 
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Figure 6: Explosion of a disc - shaped solid when only the stretched beams 
(e > 0) are allowed to break. Snapshots of the evolving system are pre- 
sented when the initial compression wave reaches the boundary of the solid 
{t = 0.0001s), the constructive interference of the incoming and outgoing 
elongation waves breaks the boundary layer (t = 0.001s) and the final break- 
ing scenario {t = 0.004s). The color scale shows the color code for the beams. 
Black indicates that the strain is close to the stretching threshold t^. 
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Figure 7: Propagation of the elastic wave in the disc - shaped sohd. (a) 
after t = lO^^sec of the initial hit, (b) after t = 3 ■ 10~^sec the compression 
wave is approaching the boundary and (c) after t = 5 ■ 10"*^ the constructive 
interference of the incoming and the outgoing elongation waves. 
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Figure 8: Explosion of a diSe.^iW.^^^-j^r^^j^ 
compression and elongation. T^^'^al 
in Fig.|. 



fallowing the beam breaking under 
rots are taken at the same times as 
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m (g) 



Figure 9: Fragment mass histogram allowing the beam breaking under com- 
pression and elongation. The breaking thresholds were chosen = 3% and 
te = 4". (a) shows on a logarithmic plot that F{m) does not exhibit power 
law behavior and an exponential function seems to be a reasonable fit (bold 
line), (b) demonstrates the quality of the fit on a semi - logarithmic plot. 
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Figure 12: The 
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Figure 13: The -fragment mass histograrniflvaryii 
te = 1° — 7°. "the stretching threshold is fixcPn 
the single pdl^gcSis are ignored. For increasing ttie' 
which is deteritii aed by the, fixed stretchiifg mj^e. 
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Figure 15: Fragmentation of a concrete block composed of elastic grains. One 
grain at the lower middle part of the block is given a high velocity directed 
inside the block. Here the velocity was 400m/s. The size of the block was 
chosen to be 40cm x 20cm. Snapshots of the evolving system are presented 
at t = 0.0004s, t = 0.0015s and t = 0.003s. 
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